###############################################
# Function to calculate SDI2 and UDI2
# by Heather J. Gunn
# Created 10/08/19
# Last modified 10/14/19
###############################################

###############################################
# Function
###############################################

SDI2.UDI2 = function(p, loading1, loading2, intercept1, intercept2, stdev2, fmean2, fsd2,
                     nodewidth = 0.01, lowerLV = -5, upperLV = 5) {
  # p = number of indicators
  # loading1 = factor loading(s) for group 1
  # loading2 = factor loading(s) for group 2
  # intercept1 = measurement intercept(s) for group 1
  # intercept2 = measurement intercept(s) for group 2
  # stdev2 = observed standard deviation(s) of the p indicators for group 2
  # fmean2 = factor mean for group 2
  # fsd2 = factor standard deviation for group 2
  # nodewidth = space between nodes during quadrature approximation (default = .01)
  # lowerLV = lowest latent variable value evaluated (default = -5)
  # upperLV = greatest latent variable value evaluated (default = 5)
  
  
  # Define evaluation of latent variable
  LV = seq(lowerLV,upperLV,nodewidth)
  
  # Create empty matrices for future arrays, matrices, etc.
  DiffExpItem12 = matrix(NA,length(LV),p)
  pdfLV2 = matrix(NA,length(LV),1)
  SDI2numerator = matrix(NA,length(LV),p)
  UDI2numerator = matrix(NA,length(LV),p)
  SDI2 = matrix(NA,p,1)
  UDI2 = matrix(NA,p,1)
  
  # Calculate SDI2 and UDI2
  for(j in 1:p){
    for(k in 1:length(LV)){
      # Calculate difference in expected indicator scores between groups 1 and 2
      DiffExpItem12[k,j] <- (intercept1[j]-intercept2[j]) + (loading1[j]-loading2[j])*LV[k]
      # probability density function for sample estimate of group 2 latent variable distribution
      pdfLV2[k] = dnorm(LV[k], mean=fmean2, sd=fsd2)
      
      # Multiply by latent variable distribution to calculate individual data point in numerator
      SDI2numerator[k,j] = DiffExpItem12[k,j]*pdfLV2[k]*nodewidth
      UDI2numerator[k,j] = abs(SDI2numerator[k,j])
    }
    # Sum across range of latent variable using quadrature to calculate numerator & divide by denominator
    SDI2[j] <- sum(SDI2numerator[,j])/stdev2[j]
    UDI2[j] <- sum(UDI2numerator[,j])/stdev2[j]
  }
  colnames(SDI2) = "SDI2"
  colnames(UDI2) = "UDI2"
  
  return(list(SDI2=round(SDI2,4), UDI2=round(UDI2,4)))
}

###############################################
# Examples
###############################################

# User-supplied parameter estimates based on output for one indicator
# p = 1
# loading1 = .8
# loading2 = .5
# intercept1 = 0
# intercept2 = -0.4
# stdev2 = 1
# fmean2 = -.5
# fsd2 = 1
# 
# effectSize = SDI2.UDI2(p, loading1, loading2, intercept1, intercept2, stdev2, fmean2, fsd2)
# effectSize$SDI2
# effectSize$UDI2

# User-supplied parameter estimates based on output for multiple indicators
# p = 8
# loading1 = c(0.9, 0.8, 0.7, 0.6, 0.8, 0.8, 0.8, 0.8)
# loading2 = c(0.9, 0.8, 0.8, 0.7, 0.7, 0.4, 0.7, 0.4)
# intercept1 = c(0,0,0,0,0,0,0,0)
# intercept2 = c(0,0,0,0.2,-0.4,0.4,-0.6,-0.6)
# stdev2 = c(1,1.5,2,1.7,1,1.3,5,6)
# fmean2 = -0.5
# fsd2 = 1
# 
# effectSize = SDI2.UDI2(p, loading1, loading2, intercept1, intercept2, stdev2, fmean2, fsd2)
# effectSize 
